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The  influence  of  attractive  absorbate-absorbate  interactions  and  surface  diffusion  on  phase  transiuons  in  ummolecular  surface 
reactions  is  examined  by  Monte  Carlo  simulations.  Detailed  studies  in  parameter  space  are  performed  assuming  various  forms  of 
desorpuon  and  reacuon  rate  expressions.  Introduction  of  rate  anomalies  ran  lead  either  to  a  kinetic  or  to  an  equilibrium  phase 
transition.  depending  on  the  reaction  luneucs  and  the  rate-determining  proce  ss.  It  is  demonstrated  that  migration  affects  the 
reaction  rate,  nudeation.  and  phase  transitions,  and,  therefore,  detailed  study  of  them  at  a  molecular  level  is  necessary  to  describe 
heterogeneous  reactions.  Hysteresis  in  reactant  coverage  and  reaction  rate  are  studied  and  large  fluctuations  are  found  to  dominate 
the  system  dynamics  near  cusp  points.  The  influence  of  defects  which  act  as  nudeation  centers  on  metastability  is  also  investigated. 
The  Monte  Carlo  data  are  compared  with  the  predictions  of  mean-field  theory,  and  limitations  of  the  latter  are  discussed. 


1.  Introduction 

Surface  reaction  systems  frequently  exhibit 
complicated  dynamical  behavior  (1-7];  reaction 
rate  multiplicities,  rate  oscillations,  phase  transi¬ 
tions,  etc.  These  non-equilibrium  systems  are  out¬ 
side  of  classical  statistical  mechanics,  and  no 
rigorous  theories  exist  to  describe  them  (6.7],  A 
classical  approach  for  the  study  of  catalytic 
processes  is  to  use  the  mean-field  approximation 
(MFA)  by  writing  a  system  of  ordinary  differen¬ 
tial  equations  in  which  a  uniform  concentration  of 
species  over  the  surface  is  assumed  (1,2,8-11). 
However,  the  MFA  cannot  handle  correctly  the 
formation  of  islands  formed  on  catalyst  surfaces 
(12,13).  In  addition,  surface  migration  is  inher¬ 
ently  excluded  from  the  model  since  atoms  are 
assumed  uniformly  distributed  and  no  concentra¬ 
tion  gradients  exist.  Generalization  of  Fick’s  Sec¬ 
ond  Law  to  account  for  adsorption,  desorption. 


This  research  partly  supported  by  NSF  under  gram  No. 
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and  reaction  with  appropriate  boundary  condi¬ 
tions  cannot  predict  uphill  diffusion  observed  in 
these  situations  (14).  As  a  consequence.  MF  mod¬ 
els.  aside  from  being  in  many  cases  analytically 
intractable  because  of  their  complexity,  cannot 
explain  the  observations  of  enhanced  operation 
under  periodic  conditions  or  the  causes  of  syn¬ 
chronization  of  multiple  oscillators  on  the  catalyst 
surface  to  produce  global  oscillations  (4,15,16). 

Monte  Carlo  (MC)  simulations  (17)  can  treat 
processes  occurring  on.  the  surface  taking  into 
account  the  spatial  inhomogeneities  and  local 
fluctuations  of  concentrations.  Recently,  Ziff, 
Guiari,  and  Barshad  proposed  a  MC  model 
(hereafter  called  ZGB  model)  to  study  a  bimolecu- 
lar  reaction  between  a  diatomic  molecule  A2  (rep¬ 
resenting  O,)  and  a  non-dissociative  molecule  B 
(representing  CO)  [6).  In  addition,  a  bimoiecular 
reaction  between  an  atom  A  and  an  atom  B 
(hereafter  called  AB  model)  has  been  studied  (18). 
Subsequently,  extensions  of  the  ZGB  model  and 
comparison  with  the  MF  theory  have  been  de¬ 
scribed  by  many  groups  (5,19-24).  Even  though 
these  models  simplify  the  experimental  situations. 
they  provide  insight  into  irreversible  phase  transi¬ 
tions  for  conditions  far  from  equilibrium. 
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In  this  work,  the  MC  method  is  used  to  ex¬ 
amine  unimolecular  catalyzed  reactions  because 
more  thorough  investigation  in  parameter  space  is 
possible  compared  with  bimolecular  reactions 
[10,25].  Furthermore,  we  suggest  that  surface  uni¬ 
molecular  reactions  involve  many  features  which 
are  present  in  real  bimolecular  systems.  For  exam¬ 
ple,  if  0  is  the  surface  coverage,  the  (1  -6)  term 
used  here  can  be  a  simplified  representation  of  the 
coverage  of  a  second  species  with  competitive 
adsorption,  i.e.  a  site  not  occupied  by  a  reactant 
can  be  vacant  or  blocked  by  a  second  species. 

Recently,  we  used  the  MC  method  to  investi¬ 
gate  the  problem  of  kinetic  oscillations  for  a  uni¬ 
molecular  reaction  in  a  continuous  stirred  tank 
reactor  [7].  In  this  paper,  the  phenomenon  cf 
metastability  for  conditions  close  to  equilibrium 
and  the  phase  transition  caused  by  reaction  vacant 
sites  are  studied  for  a  unimolecular  reaction  in  a 
constant  pressure  reactor.  Since  surface  diffusion 
cannot  be  modeled  by  a  MF  model  but  it  is  the 
dominant  process  in  many  experiments,  its  effect 
on  phase  transitions  and  nucleation  is  considered. 

It  is  observed  experimentally  that  all  surfaces 
contain  defect  sites  such  as  dislocations,  steps, 
vacancies,  etc.  Consequently,  the  binding  energy 
of  atoms  to  the  substrate  is  not  necessarily  con¬ 
stant  over  all  the  surface  [26].  At  present,  no  clear 
picture  of  the  influence  of  defects  or  impurities  on 
global  oscillations  exists  [18.27],  In  the  studv  of  a 
continuous-flow  reactor,  we  have  shown  that  de¬ 
fects  might  strongly  affect  the  existence  of  kinetic 
oscillations  [7]  through  the  influence  on  phase 
transitions.  Therefore,  their  role  on  the  dynamics 
of  heterogeneous  reactions,  which  is  absent  from 
the  MF  model,  must  be  appropriately  examined. 
In  this  paper,  we  calculate  the  isotherms  on  imper¬ 
fect  surfaces  and  investigate  how  defects  affect 
nucleation.  metastabilitv.  and  the  location  of  the 
transition  point. 

2.  Models 

2.  J.  Mean  field  model 

In  »he  MFA.  the  adatoms  are  assumed  uni¬ 
formly  distributed  on  the  surface.  In  an  isother¬ 


mal.  mixed  reactor,  the  overall  adsorption- 
desorption- reaction  process  is  described  by  [10] 


where  r3,  rd,  rR  are  the  adsorption,  desorption, 
and  reaction  rates,  and  0  is  the  surface  coverage, 
i.e.  the  number  of  atoms  adsorbed  on  the  surface 
divided  by  the  number  of  surface  sites.  The  gas 
phase  reactant  partial  pressure  P  is  assumed  to  be 
uniform  over  the  surface  and  constant  in  time  (a 
differential  reactor). 

In  the  presence  of  absorbate  attractive  interac¬ 
tions  of  strength  w<  w  >  0),  the  change  of  reactant 
coverage  with  time  is  given  in  MFA  by 

~  =  mak3P(\ -0)m‘ 

-  m3kd0  extp{~zi0w/kT)8m‘ 

-kR6(\-8)m\  (2) 

where  is  either  one  or  two  for  first-  or  second- 
order  kinetics  in  desorption,  r$  is  the  surface 
coordination  number,  and  mR  is  the  number  of 
vacant  sites  required  for  reaction.  ma  is  the 

desorption  rate  in  the  low-coverage  limit.  kR  is 
the  reaction  rate  constant,  and  k3  is  the  adsorp¬ 
tion  rate  constant  over  the  density  of  surface  sites, 
nn.  For  ideal  gas  kinetics  k3  is 

k  - (3) 

(2irmkT)l/2 

where  s  is  the  sticking  coefficient,  m  is  the  mass 
of  an  atom,  k  is  the  Boltzmann  constant,  and  T  is 
the  temperature. 

According  to  the  MFA.  each  atom  sees  a  mean 
field,  an  average  number  of  nearest  neighbors 
equal  to  z$.  A  factor  of  the  form  exp<  —  zf)w/ 
kT)  in  the  reaction  rate  of  eq.  (2)  appears  to  be 
equivalent  to  the  same  factor  in  the  desorption 
rate  which  is  examined  in  sections  3-7.  However, 
this  term  would  couple  differently  with  the  gas 
pressure  equation  in  the  modeling  of  a  flow  reac¬ 
tor  [7.10].  and  this  case  does  not  exhibit  oscilla¬ 
tions  [10].  Therefore,  in  the  calculations  of  a  con¬ 
stant  pressure  reactor  presented  here,  the  reaction 
activation  energy  is  taken  to  be  coverage  indepen- 
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dent.  The  nonlinear  terms  in  desorption  and  reac¬ 
tion  rates  of  eq.  (2)  express  the  reduction  of  ef¬ 
ficiency  in  adatom  removal  caused  either  by  at¬ 
tractive  interactions  or  vacant  sites  required  for 
example  by  a  large  organic  molecule  to  the  decom¬ 
pose  into  products.  In  reality,  these  nonlineanties 
can  act  independently  or  simultaneously.  In  most 
of  our  calculations,  the  reaction  rate  is  taken  to  be 
proportional  to  the  reactant  coverage  (rR  =  kK0). 
During  growth  of  clusters,  the  desorption  rate  can 
be  reduced  by  a  few  orders  of  magnitude  (depend¬ 
ing  on  w/kT)  compared  with  the  reaction  rate 
which  increases  with  coverage.  Hence,  the  impor¬ 
tance  of  the  two  processes  in  removal  of  atoms 
during  the  transient  path  of  the  system  can  vary 
considerably. 

If  =  l,  mR  =  0.  and  w/kT  =  0.  these  kinet¬ 
ics  are  called  Langmuir  kinetics  and  the  MF  model 
is  exact.  However,  inclusion  of  absorbate-ab- 
sorbate  interactions  or  vacant  site  requirements 
for  either  adsorption  or  reaction  can  cause  devia¬ 
tions  from  Langmuir  kinetics,  and  in  these  situa¬ 
tions  the  MF  model  is  not  exact.  In  the  following 
sections,  the  nonlinear  terms  in  desorption  and 
reaction  rates  are  examined  and  comparison  of  the 
MF  model  with  the  MC  calculations  is  presented. 

2.2.  Monte  Carlo  model 

The  surface  is  modeled  as  a  two-dimensional 
square  lattice  with  periodic  boundary  conditions. 
Each  time  a  site  is  randomly  chosen.  If  this  site  is 
empty,  adsorption  is  attempted  with  probability 
pa.  If  the  site  is  occupied,  the  atom  can  either 
desorb,  react,  migrate  with  probabilities  pd,  pK , 
and  pm  respectively,  or  not  move.  To  choose 
between  these  events,  a  segment  proportional  to 
the  probability  of  each  process  is  taken  and  a 
random  number  is  used  to  select  which  event  will 
be  executed.  To  perform  the  stochastic  simulation, 
one  has  to  describe  the  local  coverage  dependence 
of  the  probabilities  of  the  vanous  processes. 

The  basic  steps  in  the  heterogeneous  catalysis 
of  unimolecular  reactions  consist  of  adsorption  of 
the  reactant  on  the  surface,  desorption  of  the 
reactant  and  of  the  product(s)  of  the  reaction,  and 
the  reaction  step  (a  Langmuir- Hinshelwood  pro¬ 


cess).  Adsorption  of  the  reactant  A  is  assumed  to 
occur  by  chemisorption  to  a  vacant  site  S 

A  +  S  -  AS,  (4) 

where  AS  denotes  the  adsorbed  species.  Each  time 
an  atom  finds  an  empty  site  it  adsorbs  with  prob¬ 
ability  p3  which  is  proportional  to  its  gas  partial 
pressure  P ,  pa  a  kiP. 

Desorption  of  reactant  A  can  occur  back  to  the 
gas  phase  with  probability  proportional  to  the 
desorption  rate 

AS  — »  A  +  S.  (5) 

liberating  one  vacant  site.  With  attractive  interac¬ 
tions  between  the  atoms,  the  activation  energy  for 
desorption  is  coverage  dependent  and  the  prob¬ 
ability  for  desorption  is 

Pd.l<xkioexr{-'Lwjc/kT)'  (6) 

where  w  is  the  magnitude  of  interaction  between 
an  atom  on  site  i  and  its  y  th  neighbor  and  cy  is 
the  value  of  the  occupation  function  at  site  j 
which  is  either  zero  if  the  site  is  empty  or  one  if 
the  site  is  occupied.  The  reactant  coverage,  0.  is 
the  the  average  value  of  the  occupation  function 
over  all  the  surface 


(7) 


where  is  the  number  of  surface  sites.  If  the 
atoms  are  uniformly  distributed,  i.e.  the  value  of 
the  occupation  function  is  constant  over  all  the 
surface,  the  desorption  probability  described  by 
eq.  (6)  coincides  with  that  given  by  the  MF  model 
in  eq.  (2).  Here,  we  shall  only  consider  first  nearest 
neighbor  interactions  of  strength  w. 

The  rate  of  desorption.  rd.  is  calculated  in 
monolayers  per  unit  time  t 


,  d(Ad/iVJ 

rj  -  *do  d, 


(8) 


where  •Yi  the  number  of  atoms  desorbed  at  time 

t. 

Reaction  of  adsorbed  species  A  to  produce 
product  B  occurs  with  probability  proportional  to 
the  reaction  rate  constant  A:R 


AS  -»  B  +  S. 


(9) 
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The  product(s)  B  formed  is  assumed  to  desorb 
instantly  and  not  interact  further  with  the  system. 
If  two  vacant  sites  (mR  =  2  in  the  MFA)  are 
required  for  deposition  of  the  reaction  products  of 
a  polyatomic  molecule  A,  i.e. 

AS  +  2S-B  +  3S  (10) 

then  the  probability  for  reaction  of  an  atom  on 
site  t  is 


Pu..t  a  k n >. 1 1  £'j)0  £’*)> 


(ID 


where  r  ck  are  the  values  of  the  occupation 
function  at  nearest  neighbor  sites  j  and  k.  Sites  j 
and  k  are  randomly  chosen.  Thus  the  average 
probability  for  reaction  is  proportional  to  k  R, 
k r/2,  and  kR/ 6  if  there  are  4.  3,  and  2  empty 
nearest  neighbors  respectively.  Similarly,  for  the 
case  of  one  vacant  site  requirement,  only  one 
nearest  neighbor  site  is  randomly  checked  to  de¬ 
termine  the  value  of  its  occupation  function. 

In  the  MC  simulation  the  reaction  rate  is 
calculated  from  the  number  of  atoms  produced. 
•VB.  «.e. 


,  d  (Nb/NJ 

rR~*R  df 


(12) 


In  our  model,  surface  diffusion  to  one  of  the 
four  nearest-neighbor  surface  sites  is  also  allowed, 
invoking  the  principle  of  microscopic  reversibility 
or  detailed  balance  (28.29]  for  each  diffusion  jump. 
In  general,  the  dynamics  of  thermally  excited  phe¬ 
nomena  cannot  be  described  by  Kawasaki  dy¬ 
namics  or  Metropolis  walk  in  which  the  hopping 
probability  depends  only  on  the  energy  difference 
between  the  final  and  the  initial  state  of  the  sys¬ 
tem  (29].  Actually,  the  migration  probability  would 
depend  on  the  energy  barrier  for  an  atom  to  hop 
(29].  Following  Gilmer  and  Bennema  >28],  he 
migration  rate.  pm  i,  of  an  atom  on  site  t  jumping 
to  site  j  is 


Pm .  i  ^  PtS .  i 


(1  -<■,)• 


(13) 


where  xs  is  the  mean  displacement  of  an  isolated 
atom  during  its  time  life  on  the  substrate  and 
xja  is  the  average  number  of  jumps  executed  by 
an  atom  on  the  surface  before  it  desorbs.  In  this 
model,  the  energy  barrier  depends  only  on  the 


local  environment  of  site  /  (28,30]  throcch  pd ,,  eq. 
(6). 

Most  of  the  simulations  have  been  performed 
for  a  30  X  30  surface  since  size  effect  analysis 
indicated  no  appreciable  effect  if  the  lattice  size 
was  greater  or  equal  than  the  one  used.  Simula¬ 
tions  were  run  on  DEC.  Apollo.  Sparc  worksta¬ 
tions,  and  the  Cray  supercomputers  of  the  Min¬ 
nesota  Supercomputer  Institute.  Close  to  a  transi¬ 
tion  point,  a  typical  run  of  108  attempts  with  a 
program  including  adsorption,  desorption,  reac¬ 
tion,  and  diffusion  requires  approximately  30-45 
minutes  of  Cray  CPU  time  (without  vectorization) 
and  8-10  CPU  hours  on  an  Apollo  DN4500. 

2.3.  Real  lime  and  Monte  Carlo  trials 


Every  attempted  event  corresponds  to  a  Monte 
Carlo  trial  (MCT).  Null  trials  are  also  allowed  in 
the  model.  One  Monte  Carlo  step  or  sweep  (MCS) 
consists  of  Ns  MCT.  i.e.  it  is  the  time  required  to 
sample  each  surface  site  once  on  average.  MCS  is 
the  usual  reported  units  of  time  in  the  literature 
(6.18].  To  compare  the  MC  data  with  experimental 
results,  the  relationship  between  real  time  and 
MCT  is  required.  This  is  accomplished  using  the 
adsorption  process  [31]  since  no  interactions  have 
to  be  calculated  over  all  the  surface  at  each  mo¬ 
ment.  If  all  vacant  sites  for  adsorption  were  con¬ 
sidered  to  be  equivalent,  then  the  average  prob¬ 
ability  per  unit  time  for  adsorption  on  the  surface 
is 

i  * 

Pa.av=  V  -C,)-ktP(l-9).  (14) 

'  5  ,-l 


Each  time  a  successful  adsorption  event  occurs, 
the  time  elapsed  A/  corresponds  to  the  average 
adsorption  lime  per  site.  i.e. 


(1 


-  0  )  - 


(15) 


where  =  1  /ka  P  is  the  adsorption  time  on  a 
clean  surface.  Since  the  time  elapsed  between  ad¬ 
sorption  events  is  random  [32].  the  time  step  is 
then 


A  r  =  - 


(1  -9)NS 


ln(  £ ), 


(16) 


where  £  is  a  random  number  between  0  and  I. 
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The  total  elapsed  time  for  the  simulation  is  com¬ 
puted  from 


age  as  the  reactant  pressure  rises,  fig.  la.  Fig.  lb 
shows  a  typical  plot  of  the  reactant  coverage  versus 
time  obtained  by  the  MC  method  starting  either 


where  i  runs  over  all  time  steps.  When  the  MF 
model  is  exact,  calculation  of  the  time  elapsed  to 
reach  a  certain  coverage  indicates  agreement  be¬ 
tween  the  two  models  [33]. 


3.  Phase  transitions  and  metastabilitv  without 


In  two-dimensions,  the  MF  theory  is  qualita¬ 
tively  correct  in  predicting  a  coverage  multiplicity 
(three  solutions)  for  w/kT>  1  in  the  limit  of  no 
reaction,  as  shown  in  fig.  la  for  w/kT  =  2.  How¬ 
ever,  this  critical  value  of  temperature  (above 
which  multiplicity  is  observed)  is  overestimated 
and  the  exact  value  at  the  transition  point  is 
known  to  be  1.76  [17.341.  In  the  MF  theory  two  of 
the  three  solutions  are  stable  (S)  and.  once  the 
system  is  in  one  of  them,  it  will  remain  in  it  for  an 
infinite  time.  An  improvement  over  the  MF  theory 
is  accomplished  by  considering  pairs  of  sites  on 
the  surface  with  first  nearest  neighbor  interactions 
between  adatoms  instead  of  uniformly  distributed 
atoms.  The  isotherm  obtained  [35],  shown  also  in 
fig.  la.  is  called  Fowler-Guggenheim  (FG)  and  is 
given  by 


1  +  26 
26 


0.0 
0.01 


*  '\ 

e  =  i  \ 


=  o  \ 


-0„  =  0.85  , 


eD  =^o,5 


exp(  -zjv/kT). 


fee.  = 


where  K  is  the  adsorption-desorption  equilibrium 
constant.  K  =  k d/kM),  and  b  is  defined  as 

b=  {1  +40(1  -6)[exp[w/kT)  -  1)]  }'/2.  (19) 

The  characteristics  of  the  FG  isotherm  concerning 
multiplicity  and  stability  of  the  solutions  are  the 
same  with  those  of  the  MF  isotherm.  However,  the 
agreement  with  the  MC  data  is  considerably  im¬ 
proved.  The  multiplicity  regime  is  large  for  the 
MF  model  and  smail  for  the  MC  method  with  the 
FG  in  between. 

The  MC  method  indicates  an  abrupt  (first- 
order)  transition  [17]  from  a  low  to  a  high  cover- 


Fig.  1.  Panel  (a)  shows  the  FG  isotherm,  the  MF  isotherm,  and 
the  MC  data,  along  with  their  corresponding  tie  lines.  Metasta¬ 
ble  states  are  indicated  in  panel  (b)  which  shows  typical 
reactant  coverage  time  senes,  starting  from  a  clean  (0o  —  O. 
KP  -  0.0195)  or  a  totally  covered  surface  ( 9q  =  1.  KP  **  0.0172). 
and  snapshots  dunng  nuclei  growth.  Peaks  in  the  tune  series 
which  are  observed  before  the  transition  correspond  to  growth 
10O  =  U>  or  shnnktng  (ft,  =  1 )  of  nuclei  which  subsequently 
shrink  or  grow.  Panel  (c)  shows  the  reactant  coverage  as  a 
function  of  time  for  various  values  of  initial  coverage,  8a  and 
KP  ^  0  0186.  Initially,  the  atoms  are  randomly  deposited  on 
the  surface  except  for  the  case  t)„  =  0.5  where  the  atoms  are 
either  randomly  distributed  ( 9  decreases  with  lime)  or  in  a 
stripe  (6  rises  with  timei.  The  parameters  used  are  w/kT  -  2. 

fcu  =  0.  =  1. 
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from  a  clean  surface  ( 60  =  0)  or  a  totally  filled 
surface  ( 0a  -  1 )  where  0Q  is  the  initial  reactant 
coverage.  This  figure  indicates  that  without  reac¬ 
tion  metastable  states  (31]  are  found  in  which  the 
system  is  trapped  for  a  relatively  long  time  before 
it  can  overcome  the  activation  barrier  to  reach  the 
stable  equilibrium  state.  As  the  reactant  pressure 
increases,  the  transition  of  the  system  from  the 
metastable  state  (low  coverage)  to  the  stable  state 
(high  coverage)  occurs  in  shorter  time.  At  the 
equilibrium  pressure,  the  lifetime  of  the  metasta¬ 
ble  state  is  infinite  (17].  Starting  from  a  high 
partial  pressure  and  decreasing  the  reactant  pres¬ 
sure.  hysteresis  loop  is  formed,  fig.  la.  Hence 
coverage  multiplicity  is  found  by  MC  simulations 
for  an  isothermal  system. 

The  discontinuous  transition  from  a  low  to  a 
high  coverage  is  caused  by  cluster  formation  and 
growth  due  to  attractive  absorbate-absorbate  in¬ 
teractions.  Small  clusters  are  formed  by  local  vari¬ 
ation  in  the  adsorption  of  the  incoming  species 
tv  a  they  lower  the  system  free  energy.  These 
islands  grow  and  shrink  until  a  critical  cluster  is 
formed.  Then  the  probability  of  this  cluster  to 
grow  is  on  the  average  50%.  A  supercritical  cluster 
(36]  is  formed  when  one  more  atom  is  added  to 
the  critical  one.  As  the  reactant  pressure  increases, 
the  size  of  the  biggest  cluster  increases  and  when 
one  or  more  supercritical  clusters  are  formed  on 
the  surface,  they  grow  on  the  average  and  merge 
to  big  clusters  covering  almost  the  entire  surface. 
Thus  a  large  enough  perturbation  due  to  cluster 
formation  can  drive  the  system  from  the  metasta¬ 
ble  to  its  true  equilibrium  state.  Snapshots  of  the 
surface  configuration  when  the  clusters  begin  to 
grow  and  merge  are  shown  in  the  inset  of  fig.  lb. 
Since  the  transition  found  is  an  equilibrium  one, 
this  transition  is  hereafter  called  equilibrium  tran¬ 
sition.  The  effect  of  initial  conditions,  migration, 
reaction  rate  expression,  and  defects  on  the  meta¬ 
stability  are  now  discussed. 

3.1.  Dependence  on  initial  conditions 

Outside  the  multiplicity  regime,  a  unique  solu¬ 
tion  of  global  minimum  free  energy  is  obtained 
which  is  independent  of  the  initial  surface  cover¬ 
age.  However,  in  the  reactant  pressure  regime 


where  two  states  with  minima  in  the  energy  exist, 
the  state  reached  by  the  system  during  the  compu¬ 
tational  time  depends  on  the  initial  conditions. 
According  to  the  MF  theory  for  coverage  below 
the  unsteady  branch  (U),  the  low-coverage  solu¬ 
tion  is  obtained  and  for  initial  coverage  above  the 
unsteady  line,  the  high-coverage  solution  is  found. 
However,  in  the  MC  simulations,  since  the  distri¬ 
bution  of  atoms  is  not  uniform  due  to  attractive 
interactions,  the  final  state  approached  by  the 
system  during  the  computational  time  depends 
not  only  on  the  value  of  the  initial  coverage,  0o, 
but  also  on  the  distribution  of  atoms,  fig.  lc.  In 
the  MF  or  FG  approximation,  the  tie  lines  shown 
in  fig.  la.  can  be  obtained  by  the  Maxwell  equal 
area  rule  (10].  However,  in  practice  the  unsteady 
steady  branch  is  not  known.  At  the  equilibrium 
pressure,  both  states  are  stable,  i.e.  the  depth  of 
the  minima  in  the  free  energy  of  the  system  are 
equal. 

Because  both  states  are  equally  accessible,  one 
can  start  with  two  equal  stripes,  one  filled  and  one 
empty  [6],  and  monitor  the  lime  evolution  of  the 
system  until  the  surface  is  either  filled  or  almost 
empty.  Performing  30  individual  runs,  the  equi¬ 
librium  pressure  for  which  the  probability  of  going 
up  or  down  is  one  half  was  found  to  be  0.0183  ± 
0.0001  for  w/kT  =  2.  For  the  values  of  reactant 
pressure  0.0182  and  0.0184  the  probabilities  of 
going  up  are  30%  and  70%  respectively.  The  same 
simulations  were  repeated  by  locating  initially  a 
square  cluster  at  the  center  of  the  surface.  No 
effect  of  the  shape  of  the  cluster  on  the  value  of 
the  equilibrium  pressure  has  been  found.  The  MC 
points  connected  by  the  dashed  line  in  fig.  la 
correspond  to  one-half  probability  of  going  up  or 
down  starting  with  half  the  surface  covered  in  a 
stripe. 


4.  Effect  of  diffusion  on  metastability 

Surface  migration  is  another  driving  force  for 
nucleus  formation  and  growth  because  single 
atoms,  executing  many  jumps  on  the  surface  be¬ 
fore  they  desorb  (large  x./a)  have  a  large  prob¬ 
ability  during  their  random  walk  to  be  trapped  in 
a  cluster.  Atoms  in  the  interior  of  a  cluster  are  not 
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Fig.  2.  Effect  of  surface  diffusion  on  the  dynamics  of  the 
equilibrium  phase  transition  for  a  slow  reaction  (m„  =  2. 
kK/kM  =  10' 31  as  a  function  of  the  average  number  of  jumps. 
xs/a.  Panel  (a)  corresponds  to  ft,  =  0.  KP  =  0.02  and  panel  (b) 
to  ft,  =  l.  KP  =  0  0168.  Surface  diffusion  draws  the  system 
faster  from  a  low  to  a  high  coverage  in  the  presence  of 
attractive  interactions. 


allowed  to  jump  to  neighbor  positions,  unless  there 
is  a  hole  adjacent  to  them.  Furthermore,  atoms  at 
the  border  of  clusters  must  break  all  bonds  with 
their  nearest  neighbors  to  overcome  the  migration 
barrier  and  diffuse.  The  more  compact  is  the 
cluster,  the  fewer  atoms  can  escape  from  it.  Thus 
as  time  proceeds,  larger  clusters  are  formed  on  the 
surface.  Therefore,  migration  facilitates  nucleation 
and  brings  the  system  faster  to  its  stable  state,  as 
shown  in  fig.  2a.  Each  curve  is  an  average  over 
five  individual  runs. 

In  contrast  to  the  transition  from  a  low  to  a 
high  coverage,  if  one  starts  with  an  initially  covered 
surface,  at  a  sufficiently  low  reactant  pressure, 
desorption  will  create  holes  and  disintegrate  the 
surface  into  small  clusters.  Hence,  desorption 
competes  with  diffusion  which  brings  atoms  to  the 
clusters  and  the  simulations  indicate  that  diffusion 


results  in  a  slowing  down  of  the  transition  from  a 
high  to  a  low  coverage,  as  shown  in  fig.  2b. 

Therefore,  surface  migration  speeds  up  the 
transition  from  a  low  to  a  high  coverage  but  slows 
down  the  transition  from  a  high  to  a  low  coverage. 

5.  Effect  of  reaction  rate  expression  on  the  phase 
transition 

In  section  3.  the  phase  transition  was  examined 
for  the  adsorption-desorption  process.  In  this  sec¬ 
tion.  the  effect  of  reaction  on  the  phase  transition 
driven  by  attractive  interactions  is  discussed.  The 
reaction  rate  is  given  by  rR  =  kR0(  1  -  0)m*.  where 
mR  is  the  number  of  vacant  sites  required  for 
decomposition  of  a  polyatomic  molecule  A.  The 
factor  ( 1  -  0 )  can  be  also  considered  as  the  cover¬ 
age  of  a  second  species  in  a  bimolecular  reaction. 
In  the  following,  the  cases  of  zero,  one,  and  two 
vacant  sites  are  investigated. 

5.1.  Langmuir  reaction  (mR-  0) 

In  the  interior  of  clusters,  the  probability  of  an 
atom  to  desorb  is  very  small  due  to  the  large 
number  of  bonds  it  has  to  break  with  nearest 
neighbors.  Under  these  circumstances,  the  surface 
reaction  can  be  the  dominant  mechanism  which 
removes  atoms  from  the  interior  of  nuclei,  creat¬ 
ing  holes  and  inhibiting  nucleation.  If  the  reaction 
rate  constant  is  not  very  large,  the  transition  from 
a  low  to  a  high  coverage  and  reaction  rate  remains 
abrupt,  as  shown  in  fig.  3a.  However,  since  cluster¬ 
ing  is  reduced  by  the  surface  reaction,  a  higher 
reactant  pressure  is  required  for  nucleation  to 
occur  and  therefore  the  transition  is  shifted  to  a 
higher  pressure  as  kR  is  increased.  When  the 
reaction  rate  constant  is  sufficiently  large,  nuclea¬ 
tion  and  subsequent  growth  of  clusters  is  blocked 
and  a  smooth  isotherm  is  obtained,  as  shown  in 
fig.  3  a.  Starting  from  a  totally  covered  surface 
(high  reactant  pressure)  and  decreasing  the  pres¬ 
sure.  hysteresis  loop  is  formed,  fig.  3a.  and  there¬ 
fore.  reactant  concentration  multiplicity  and  reac¬ 
tion  rate  multiplicity  are  observed.  In  fig.  3.  verti¬ 
cal  lines  depict  a  discontinuous  transition  and 
arrows  show  the  direction  of  transition  (from  a 
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low  to  a  high  coverage  or  vice  versa).  Hysteresis 
loop  is  indicated,  whenever  it  exists,  by  a  pair  of 
lines  of  the  same  symbol. 

The  reduction  of  the  metastable  regime  as  kR 
increases  is  shown  in  fig.  3c  in  which  the  pressure 
Pt  at  the  two  turning  points  of  the  S-shaped 
isotherms  is  plotted  as  a  function  of  reaction  rate 
constant,  kR.  A  discontinuous  transition  is  also 
observed  for  the  average  number  of  nearest 
neighbors,  the  percolation  order  parameter,  and 


the  average  size  of  clusters  formed  on  ne  surface 
[33], 

The  reaction  rate  (not  shown)  calculated  in  the 
MC  simulations  is  found  to  be  proportional  to  the 
reactant  coverage,  i.e.  rR//cR  =  8.  within  statistical 
errors  as  is  predicted  by  the  MF  theory.  Accord¬ 
ing  to  the  MFA,  the  desorption  rate,  rd/kM  = 
8  exp (-zs0w/kT),  plotted  as  a  function  of  cover¬ 
age  should  be  independent  of  kR.  However,  devia¬ 
tions  from  the  MFA  are  observed  in  the  desorp- 
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Fig.  3.  Effect  of  reaction  rate  expression  and  reaction  rate  constant.  kR,  on  metastabilitv  for  »  cT  -  2.  Panels  (aHc)  correspond  to 
mR  =  0.  (d)-(0  1°  m r  =  1-  and  (g)-(i)  to  mR  =  2.  Lines  with  the  same  symbol  in  panels  la).  (d).  (g),  and  (h)  show  an  hysteresis  loop 
in  coverage  or  reaction  rale.  Panels  (a),  (d).  and  (g)  show  isotherms  for  the  three  models.  Panel  lb)  shows  deviations  from  the  MFA  in 
desorption  rate  for  the  mR  =  0  model  whereas  panel  (e)  shows  deviations  from  the  MFA  in  reaction  rate  for  the  mR  -  1  modei  for 
various  values  of  kR/ki0.  Panel  (h)  shows  hysteresis  loops  in  reaction  rate  for  the  mR  “  2  model.  No  single-valued  isotherms  or  rates 
are  obtained  for  this  model.  Panels  (c).  (F),  and  d)  show  the  muluplicitv  regime  as  a  function  of  kR.  No  multiple  solutions  can  be 
found  for  the  mR  =»  0  and  mR  =  1  models  above  a  certain  value  of  kR.  In  contrast  the  multiplicity  regime  for  the  mR  =  2  model  does 

not  exhibit  a  cusp  point. 
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non  rate  which  is  not  treated  correctly  because  the 
.V1F  theory  does  not  take  into  account  the  forma¬ 
tion  of  clusters  on  the  surface,  as  shown  in  fig.  3b. 
As  the  reaction  rate  constant  increases,  the  clus¬ 
ters  become  less  compact  and  the  desorption  rate 
increases.  For  kR/kM=  1.  it  is  found  that  the 
average  number  of  nearest  neighbors  seen  by  an 
atom  does  not  deviate  much  from  the  MF  theory. 
However,  large  deviations  are  still  observed  from 
the  MF  desorption  rate  due  to  erroneous  calcula¬ 
tion  of  the  desorption  rate  based  on  a  uniform 
surface  coverage  assumption.  Examination  of 
surface  configurations  and  cluster  distribution  in¬ 
dicates  that  the  reaction  creates  holes  and  pro¬ 
duces  manv  single  atoms  in  the  regime  where  the 
maximum  of  the  desorption  rale  occurs  which 
then  increase  considerably  the  desorption  rate  over 
the  MF  estimated  value. 

5.2.  One  vacant  site  frn  R  =  I )  required  for  reaction 

In  this  case  one  vacant  site  is  required  for  the 
decomposition  of  molecule  A.  For  small  values  of 
the  reaction  rate  constant  the  transition  remains 
abrupt  but  the  transition  point  shifts  to  the  right, 
i.e.  to  a  higher  reactant  pressure,  as  shown  in  fig. 
3d.  At  higher  values  of  the  reaction  rate  constant  a 
smooth  isotherm  is  obtained.  Fig.  3f  indicates 
that,  as  the  reaction  rate  constant  increases,  the 
regime  in  which  multiple  solutions  are  observed 
becomes  narrower  and  a  cusp  point  exists  for  each 
value  ot  interaction  strength  at  which  the  two 
turning  points  of  the  isotherm  coincide. 

The  qualitative  behavior  of  the  system  is  simi¬ 
lar  to  the  Langmuir  reaction  case  ( mR  =  0)  studied 
above.  However,  the  transition  from  discontinu¬ 
ous  to  smooth  isotherm  occurs  at  higher  reaction 
rate  constant.  Since  atoms  in  the  intenor  of  com¬ 
pact  clusters  cannot  be  annihilated  by  reaction 
because  they  cannot  find  any  vacant  site  for 
decomposition,  the  reaction  requiring  one  vacant 
site  is  not  as  efficient  for  destroying  clusters  as 
the  Langmuir  reaction  examined  above.  Conse¬ 
quently.  multiple  solutions  exist  up  to  higher  reac¬ 
tion  rate  constants. 

When  one  vacant  sue  is  required,  both  the 
desorption  rate  and  the  reaction  rate  calculated  bv 
the  MC  simulations  deviate  from  the  MFA.  As  an 


example,  fig.  3e  compares  the  reaction  rate  versus 
reactant  coverage  for  the  MFA  and  the  MC  model. 
Based  on  the  MF  theory,  the  reaction  rate.  rR/kR 
=  0(1  —8).  is  a  parabolic  function  of  8.  Since 
reaction  is  considered  to  be  very  slow,  the  forma¬ 
tion  of  clusters  is  caused  mainly  by  attractive 
interactions.  Nucleation  and  growth  of  clusters 
result  in  a  low  reaction  probability  and  thus  in  a 
lower  reaction  rate  compared  with  the  case  of 
uniform  adatom  distribution  as  shown  in  fig.  3e. 
When  a  discontinuous  transition  occurs.  rR/kR  is 
small  and  corresponds  either  to  a  low  or  to  a  high 
coverage  (not  shown  in  fig.  3e  for  reasons  of 
clantvi.  When  a  smooth  isotheim  is  obtained,  fig. 
3d.  high  reaction  rate  can  be  achieved  as  shown  in 
fig.  3e.  The  maximum  occurs  at  coverage  -  0.5  in 
agreement  with  the  MF  theory.  The  higher  the 
reaction  rate  constant,  the  higher  the  randomiza¬ 
tion  of  the  surface,  and  the  more  accurate  the  MF 
theory  is. 

The  large  deviations  of  the  MC  data  from  the 
MF  can  be  improved  considerably  if  one  considers 
that  the  probability  of  an  atom  to  react  is  propor¬ 
tional  to  (l~r/4).  where  r  is  the  number  of 
nearest  neighbors  of  that  atom.  The  average  reac¬ 
tion  rate  at  steady  state  can  be  calculated  from  the 
coverage  and  the  average  number  of  nearest 
neighbors  seen  by  each  atom.  i.e.  rR/k  R  =  8  (1  — 
:/A).  The  simulations  show  very  good  agreement 
between  the  reaction  rate  calculated  that  wav  and 
the  rate  obtained  by  eq.  (12).  Thus,  the  failure  of 
the  MF  theory  is  due  to  the  wrong  averaging  over 
the  whole  surface  which  neglects  the  spatial  inho¬ 
mogeneities.  i.e.  the  cluster  formation  due  to  at¬ 
tractive  interactions. 

5.3.  Two  vacant  sites  Im R  =  2)  required  for  reaction 

In  this  case  two  vacant  sites  are  required  for  the 
decomposition  of  the  reactant  A.  Isotherms  calcu¬ 
lated  for  different  values  of  the  reaction  rate  con¬ 
stant  for  w/kT  =  2  are  shown  in  fig.  3 g.  As  the 
reaction  rate  constant  increases,  more  atoms  are 
removed  by  the  chemical  reaction  and  the  steady- 
state  value  of  the  reactant  coverage  decreases.  The 
requirement  of  two  vacant  sites  is  quite  restrictive 
so  that  atoms  have  less  probability  to  react  com¬ 
pared  with  the  mR  =  1  case.  Therefore,  no  appre- 
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Fig.  4.  Dvnamics  close  to  ihe  cusp  point  shown  in  fig.  ic.  Panels  (a)  and  (b)  show  ihe  fluctuations  in  reactant  coverage  versus  time  for 
two  different  lattice  sizes.  Panel  (c)  depicts  the  amplitude  of  fluctuations.  SO.  of  a  15  x  15  square  lattice  as  a  function  of  reaction  rate 
constant  kK.  for  various  values  of  reactant  pressure.  Very  small  fluctuations  are  found  below  a  certain  pressure.  Panel  (d)  shows  the 
amplitude  of  fluctuations  versus  kR  at  KP  -  0-022  for  three  different  lattice  sizes. 


ciable  effect  of  reaction  on  clustering  occurs,  and 
rit/^R  1S  only  slightly  decreased  as  a  result  of 
coverage  reduction.  Deviations  in  both  the  reac¬ 
tion  and  desorption  rates  are  again  observed.  Fig. 
3h  shows  hysteresis  loops  in  reaction  rate  as  a 
function  of  reactant  pressure  for  various  values  of 
the  reaction  rate  constant.  No  single-valued  solu¬ 
tion  exists  for  this  model. 

In  contrast  with  the  cases  mR  =  0  and  mR  =  1. 
the  discontinuous  transition  characteristic  of  the 
phase  transition  d'te  to  attractive-attractive  inter¬ 
actions  remains  a.  o  for  conditions  far  for  equi¬ 
librium.  In  particular,  the  multiplicity  regime  in¬ 
creases  with  the  reaction  rate  constant,  as  shown 
in  fig.  3i,  i.e.  as  we  leave  equilibrium  conditions, 
and  no  cusp  point  exists.  An  abrupt  transition 
from  a  low  coverage  (a  high  reaction  rate)  to  a 
high  coverage  (a  low  reaction  rate)  is  found  for  all 


cases  investigated.  This  phenomenon  is  caused  by 
the  cooperative  action  of  desorption  and  reaction 
with  vacant  sites  requirements  on  cluster  forma¬ 
tion  and  will  be  discussed  in  more  detail  in  section 

8. 


6.  Fluctuations  close  to  a  cusp  point 

The  transition  from  single  solution  to  multiple 
solutions  in  coverage  and  reaction  rate  is  char¬ 
acterized  by  the  existence  of  a  cusp  point  when 
=  0  or  m R  =  1 .  As  the  reaction  rate  constant 
increases,  there  is  a  value  of  the  reaction  rate 
constant  (depending  on  the  interaction  strength) 
which  separates  the  single  solution  subspace  from 
the  multiple  (three)  solution  subspace.  Using  the 
MF  theory  one  can  calculate  the  cusp  point  rigor- 
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ouslv  by  solving  lhe  system  of  the  corresponding 
equations. 

Close  to  a  cusp  point,  very  large  fluctuations 
dominate  small  areas  of  the  catalyst  surface,  as 
shown  in  fig.  4a  for  the  mR  =  0  case.  The  coverage 
fluctuates  between  low  and  high  values.  However, 
this  is  a  size  effect  as  shown  in  fig.  4b.  To  analyze 
the  dynamics  of  the  time  senes  obtained  by  the 
MC  simulations  we  use  the  correlation  integral 
method  (7,18,37],  We  have  found  [33]  that  the 
correlation  exponent  is  proportional  to  the  embed¬ 
ding  dimension  (7.18.37]  as  the  latter  is  increased, 
indicating  that  the  dynamics  are  dominated  bv 
stochastic  noise. 

To  locate  a  cusp  point  in  a  MC  simulation  for 
each  value  of  the  interaction  energy,  one  has  to 
determine  both  the  reaction  rate  constant  and  the 
gas  reactant  pressure.  First  a  size  effect  analysis  is 
used  in  analogy  with  the  equilibrium  Ising  model 
[17].  One  can  either  fix  the  reaction  rate  constant 
and  vary  the  pressure  or  vice  versa,  fig.  4c.  The 
maximum  of  the  fluctuations  in  coverage  should 
indicate  the  location  of  the  cusp  point.  Below  a 
certain  reactant  pressure  no  large  fluctuations  are 
found.  This  indicates  that  these  simulations  are 
performed  in  the  multiplicity  regime  where  a  rela¬ 
tively  large  perturbation  brings  the  system  from 
one  state  (a  low  coverage)  to  the  other  (a  high 
coverage).  Once  the  system  has  left  its  original 
state,  it  does  not  seem  to  return  to  its  initial  state, 
at  least  during  the  computational  time,  if  its  size  is 
large  enough.  Very  large  fluctuations  are  observed 
when  reaction  rate  constants  are  slightly  greater 
than  the  cusp  point  value,  i.e.  as  we  approach  the 
cusp  point  from  the  unique  solution  subspace  in 
agreement  with  the  equilibrium  adsorption-de¬ 
sorption  system  behavior  (33].  The  dependence  of 
the  location  of  the  maximum  of  fluctuations  on 
the  size  is  not  very  large  as  shown  in  fig.  4d.  For 
the  calculations  shown  in  fig.  4,  we  used  a  con¬ 
stant  number  of  MCS  of  approximately  120000 
for  all  sizes  of  the  surface. 

The  reactant  pressure  at  the  turning  points  of 
the  S-shaped  isotherm  versus  the  reaction  rate 
constant  is  also  plotted  in  fig.  3c.  From  classical 
bifurcation  theory,  one  expects  that  the  two  turn¬ 
ing  points  (simple  folds)  will  coincide  at  the  cusp 
point.  This  is  actually  observed  in  fig.  3c.  and  the 


value  of  pressure  and  reaction  rate  constant  ob¬ 
tained  by  this  method  are  in  very  good  agreement 
with  the  size  effect  analysis.  Similar  dynamic  be¬ 
havior  close  to  the  cusp  point  is  also  found  for  the 
mR  =  1  case  (not  shown). 

Our  calculations  suggest  that  as  ?  continuation 
parameter  is  varied  and  a  cusp  point  is  ap¬ 
proached  from  the  smooth  isotherm  subspace, 
large  fluctuations  can  determine  the  dynamics  of 
the  system.  This  is  quite  analogous  with  the  behav¬ 
ior  observed  by  Fichthom  et  al.  [18]  in  their 
bimolecular  AB  model  where  the  desorption  prob¬ 
ability  is  the  bifurcation  parameter.  In  their  model, 
at  zero  desorption  probability  a  discontinuous  iso¬ 
therm  is  obtained  and  as  the  desorption  probabil¬ 
ity  tends  to  zero  a  transition  from  a  smooth  iso¬ 
therm  to  a  discontinuous  one  occurs.  However, 
their  model  does  not  have  a  rigorously  defined 
“cusp  point”.  This  difference  might  explain  the 
fact  that  their  correlation  exponent  is  constant  as 
the  embedding  dimension  increases  whereas  our 
correlation  exponent  is  proportional  to  the  embed¬ 
ding  dimension. 

7.  Influence  of  defects  on  phase  transitions 

In  most  of  the  MC  studies  a  square  lattice  is 
used  and  it  was  assumed  that  all  sites  are  equiv¬ 
alent.  a  highly  idealized  situation.  Experiments 
indicate  that  atoms  are  frequently  preferentially 
adsorbed  either  at  step  sites  or  terrace  sites  [38,39], 
Recently,  the  influence  of  surface  steps  on  the 
phase  transitions  has  been  studied  by  the  MC 
method  by  examining  the  size  effect  dependence 
of  the  equilibrium  lattice  gas  model  (40}. 

We  restrict  our  studies  to  sites  with  very  strong 
binding  energy  which  are  distributed  either  regu¬ 
larly  or  randomly  on  a  flat  surface  (terrace)  in  the 
limit  of  a  very  slow  reaction  (kR/kd 0=  10-5).  It 
is  assumed  that  atoms  on  defect  sites  are  very 
strongly  bound  on  the  substrate  and  they  can 
neither  desorb  nor  react.  However,  these  atoms 
interact  with  their  first  nearest  neighbors  and  con¬ 
stitute  nucleation  centers.  Regularly  distributed 
defects  are  implemented  in  the  computer  as  a 
square  cluster  at  the  center  of  a  certain  size  surface 
area. 
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Fig.  5.  Influence  of  concentration  of  defects  on  phase  transition  and  multiplicity  regime.  Panels  (a)  and  <b)  show  hysteresis  loops  in 
reactant  coverage  and  coverage  multiplicity  respectively  as  a  function  of  defect  concentration.  Panels  (c).  (d).  (e)  show:  1.2%  defect 
concentration  on  a  27  x  27  surface.  4%  on  a  30  x  30  surface,  and  7.4%  on  a  33  x  33  surface  respectively  for  KP  ”  0.016.  at  i/t,  -  10. 


Figs.  5a  and  5b  show  the  influence  of  con¬ 
centration  of  regularly  distributed  defects  on  re¬ 
actant  coverage,  phase  transition,  and  reactant 
coverage  multiplicity.  Panels  (c)-(e)  are  snapshots 
of  the  surface  for  each  concentration  used  in  figs. 
5a  and  5b.  Defects  with  very  strong  binding  en¬ 
ergy  which  can  interact  through  first  nearest 
neighbor  attractive  forces  are  nucleation  centers 
which  facilitate  nucleation.  Consequently,  the  re¬ 
actant  pressure  for  the  transition  from  a  low  to  a 
high  coverage  decreases  and  the  multiplicity  reg¬ 
ime  is  reduced,  fig.  5b.  In  addition,  the  average 
coverage  on  the  surface  is  increased  and  the  de¬ 
sorption  rate  is  reduced  compared  with  the  perieci 
surface. 

Fig.  6a  depicts  the  isotherms  obtained  for  a  4% 
constant  concentration  of  defects  but  different 
defect  distributions  shown  in  figs.  6c -6f.  As  the 
number  of  edge  and  comer  defect  atoms  increases. 


more  defects  can  interact  with  adsorbed  atoms 
and  the  number  of  nucleation  centers  increases. 
Therefore,  the  location  of  the  phase  transition  is 
shifted  to  a  lower  reactant  pressure  and  the  meta¬ 
stability  regime  becomes  more  narrow  as  shown  in 
fig.  6b.  The  simulations  predict  that  single  defect 
sites  distributed  either  regularly  or  randomly  are 
more  efficient  in  shifting  of  the  location  of  the 
phase  transition  and  in  the  reduction  of  the  meta¬ 
stability  regime. 

8.  Phase  transition  driven  by  vacant  site 
requiremei:  '  -r  reaction 

To  understand  the  different  bifurcation  behav¬ 
ior  shown  in  figs.  3c.  f.  and  i.  the  adsorption-reac¬ 
tion  model  with  vacant  sites  in  the  absence  of 
desorption  <  Ard0  =  0)  is  first  xamined.  When  two 
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vacant  sites  are  required  for  reaction  ( m  R  =  2),  the 
model  has  an  analogy  with  the  ZGB  model  [6] 
where  two  vacant  sites  are  required  for  adsorption 
of  O,  into  empty  perimeter  sites.  Our  model  is 
also  related  to  the  “contact  process”  where  an 
interacting  particle  system  is  used  to  examine  the 
spread  of  an  infection.  In  this  model,  the  desorp¬ 
tion  rate  which  corresponds  to  our  reaction  rate  is 
proportional  to  the  number  of  empty  nearest 
neighbors  (for  a  discussion  see  [25]  and  references 
therein).  Our  mR  =  1  model  is  similar  to  the  model 
studied  in  one  dimension  by  Dickman  and 
Burschka  in  which  a  particle  reacts  if  at  least  one 
of  the  nearest  neighbors  sites  is  vacant  (25). 

Both  mR  =  1  and  mR  =  2  models  show  devia¬ 
tions  from  MFA  in  one  and  two  dimensions  [33]. 
For  brevity,  only  the  case  of  two  vacant  sites  for 
reaction  (mR  =  2)  is  presented  here  m  two  dimen¬ 
sions.  An  abrupt  phase  transition  (hereafter  called 
kinetic  transition)  from  a  low  coverage  to  a  totally 
covered  surface  where  the  reaction  rate  is  zero 


(poisoning  of  reaction)  is  observed  as  the  leactant 
pressure  increases.  This  transition  which  is  ob¬ 
served  when  two  vacant  sites  are  needed  for  re¬ 
actant  decomposition,  occurs  at  much  higher  pres¬ 
sure  than  the  equilibrium  transition  studied  above 
(sections  3-7)  and  is  not  due  to  absorbate-ab- 
sorbate  interactions  since  desorption  is  not  in¬ 
cluded  in  most  of  the  simulations  of  this  section. 
It  is  of  a  kinetic  character  because  the  overall 
process  is  irreversible.  In  the  following,  the  ad¬ 
sorption-reaction  mR  =  2  model  is  examined.  Fi¬ 
nally.  the  effect  of  Langmuir  desorption  (w/kT  — 
0 )  and  migration  on  the  transition  is  also  dis¬ 
cussed. 

The  MFA  predicts  that  for  low  pressures 
( k3P/kR  <  0.25)  three  solutions  exist,  fig.  7a.  The 
high-coverage  solution  corresponds  to  a  poisoned 
catalytic  surface  whereas  the  middle-coverage  one 
is  unstable.  The  only  reactive  attractor  corre¬ 
sponds  to  the  lower  coverage.  The  MC  data  indi¬ 
cate  a  discontinuous  transition  from  a  low-cover- 
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age  to  a  poisoned  surface  at  a  pressure  much 
lower  than  the  MF  prediction,  fig.  7a.  At  very  low 
reactant  pressure,  a  few  atoms  are  present  on  the 
surface  and.  since  they  are  almost  randomly  dis¬ 
tributed.  very  good  agreement  with  the  MF  is 
found.  As  the  reactant  pressure  increases,  local 
spatial  fluctuations  occur  which  lead  to  the  forma¬ 
tion  of  smaii  clusters.  For  rectangular  clusters, 
reaction  can  take  place  only  at  their  borders.  Sites 
in  the  interior  of  the  nuclei  are  inactive.  As  the 
reactant  pressure  increases  the  cluster  size  rises 
and  when  it  reaches  a  critical  value  the  surface  is 
rapidly  poisoned.  Starting  with  totally  covered 
surface  and  reducing  the  reactant  pressure,  no 
reactive  state  can  be  reached  without  desorption. 

The  effect  of  Langmuir  desorption  ( w/kT  =  0. 
no  diffusion)  is  shown  in  fig.  7a  for  different 
values  of  desorption  constant,  kd0/k  R.  Desorp¬ 


tion  blocks  nucieation  and  unless  the  desorption 
probability  is  much  less  than  the  reaction  one,  no 
discontinuous  transition  is  observed.  Since  desorp¬ 
tion  disintegrates  the  nuclei,  a  higher  reactant 
pressure  is  required  for  nucieation  to  occur,  i.e. 
the  transition  is  shifted  to  the  right  of  the  kinetic 
transition.  The  creation  of  holes  in  the  interior  of 
nuclei  and  the  removal  of  atoms  bound  at  the 
cluster’s  periphery  result  in  a  lower  degree  of 
clustering  on  the  catalyst  [33],  The  desorption  of 
atoms  from  the  interior  of  clusters  makes  them 
less  compact  and  more  molecules  can  find  two 
vacant  sites  to  deposit  the  reaction  products.  Con¬ 
sequently.  higher  reactivity  is  achieved,  fig.  7b. 
When  desorption  is  slow,  hysteresis  loop  is  formed 
if  we  start  with  a  totally  covered  surface.  Hence, 
desorption  is  the  main  process  against  poisoning 
of  the  catalytic  surface.  Due  to  desorption,  it  is 


Fig.  7.  Phase  transition  driven  bv  vacant  sites  requirements  for  reaction.  The  kinetic  phase  transition  for  mR  =  2  in  the  absence  of 
Jesorpuon  and  diffusion  is  shown  bv  the  dash-dotted  line  in  panels  (a)  and  tel  whereas  the  MF  model  is  shown  by  the  solid  line. 
Panels  tal  and  (b|  show  the  effect  of  Langmuir  desorption  (  w/kT  =  0)  without  diffusion  whereas  panels  (c)  and  tdf  show  the  effect  of 
surface  diffusion  un  the  limit  of  a  very  slow  desorption)  on  the  reactant  isotherm  and  reaction  rate  respectivelv. 
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possible  to  maintain  a  higher  coverage  on  the 
surface  and  to  increase  the  maximum  reaction  rate 
obtained. 

The  effect  of  surface  diffusion  on  the  kinetic 
transition  is  depicted  in  figs.  7c  and  7d.  The 
curves  at  the  left  of  the  kinetic  transition  corre¬ 
spond  to  coverage  dependent  diffusion  coefficient 
(in  the  limit  of  <cd0  -*  0,  w/kT  =  2).  When  lateral 
interactions  are  important,  the  diffusion  sets  a 
driving  force  for  providing  atoms  to  the  cluster 
resulting  in  bigger  clusters  on  the  surface.  The 
steady-state  distribution  of  clusters  indicates  that 
in  this  case  fewer  single  atoms  are  present  on  the 
surface  and  the  size  of  the  largest  cluster  increases 
with  xja  (33).  This  increase  of  coverage  and 
degree  of  clustering  on  the  catalytic  surface  leads 
to  more  easy  poisoning  of  the  catalyst,  i.e.  shifting 
of  the  transition  to  the  left  (lower  reactant  pres¬ 
sure),  fig.  7c. 

The  curves  at  the  right  of  the  kinetic  transition 
correspond  to  coverage  independent  diffusion 
coefficient  ( kd0  =  0.  w/kT~  0).  i.e.  no  interac¬ 
tions  between  the  adatoms.  This  situation  corre¬ 
sponds  to  a  “Fickian  kind”  of  diffusion.  In  this 
case  single  atoms  and  atoms  at  the  periphery  of 
clusters  have  equal  probability  for  migration.  In 
the  reactant  pressure  regime  where  clustering  is 
important  and  many  atoms  exist  at  the  periphery 
of  nuclei,  on  the  average  more  atoms  escape  from 
nuclei.  The  average  number  of  single  atoms  on  the 
surface  increases,  and  the  size  of  largest  cluster 
decreases  [33].  Therefore,  the  diffusion  process 
randomizes  the  surface.  As  the  mean  displacement 
xja  of  the  atoms  on  the  surface  increases,  previ¬ 
ously  blocked  atoms  can  react.  For  very  large 
migration  rates,  the  surface  has  been  almost 
randomized,  and  the  MFA  is  quite  accurate. 

Comparison  of  our  mR  =  2  model  with  the 
ZGB  model  indicates  that  the  feature  of  two  va¬ 
cant  sites  in  reaction  or  adsorption  causes  a  dis¬ 
continuous  transition  to  a  poisoned  surface  and 
active  steady  states  (<?*  1)  are  only  possible  in 
two  or  more  dimensions,  i.e.  in  one  dimension  the 
chain  becomes  eventually  poisoned  [33].  In  ad¬ 
dition.  a  Langmuir  desorption  process  can  destroy 
the  discontinuous  transition  and  a  Fickian  diffu¬ 
sion  randomizes  the  surface  and  shifts  the  transi¬ 
tion  to  a  higher  pressure.  These  observations  are 


in  qualitative  agreement  with  the  studies  of  the 
role  of  these  processes  on  the  ZGB  discontinuous 
phase  transition  [5,24], 

9.  Conclusions 

Detailed  Monte  Carlo  studies  have  been  used 
to  examine  the  adsorption-desorption-reaction- 
diffusion  processes  for  a  unimolecular  reaction  on 
a  catalyst  surface  using  various  forms  for  adsorp¬ 
tion,  desorption,  and  reaction  rates.  Formation  of 
clusters  and  their  subsequent  growth  has  been 
observed  in  our  two-dimensional  simulations  along 
with  their  profound  influence  on  surface  reactions. 
Phase  transitions  can  be  driven  either  by  attractive 
absorbate-absorbate  interactions  or  by  vacant 
sites  required  for  reactant  decomposition.  In  both 
cases  abrupt  first-order  transitions  have  been 
found  for  appropriate  values  of  the  parameters. 
The  existence  of  discontinuous  transitions  is  usu¬ 
ally  associated  with  hysteresis  loops  which  indi¬ 
cate  multiplicity  in  coverage  and  desorption-reac¬ 
tion  rates. 

The  MF  theory  is  qualitatively  correct  in  two 
dimensions  in  predicting  the  existence  of  multiple 
solutions  and  cusp  points.  However,  large  devia¬ 
tions  in  coverage  and  desorption-reaction  rates 
can  be  found  between  MF  theory  and  MC  simula¬ 
tions.  In  general,  the  formation  of  clusters  caused 
by  any  reason  cannot  be  treated  correctly  by  the 
MFA  because  of  its  uniform  coverage  assumption. 
Other  investigations  of  bimolecular  systems  of  re¬ 
actants  A  and  B  without  any  attractive  interac¬ 
tions  or  vacant  sites  requirements  indicate  that 
deviations  between  the  MF  model  and  MC  data 
would  be  expected  if  the  desorption  rate  constants 
of  the  reactants  are  small  so  that  the  randomiza¬ 
tion  of  atoms  on  the  surface  is  slow  [33].  Thus,  the 
nonlinearity  in  the  reaction  term,  rR  =  fcR0A0B, 
can  result  in  inaccuracy  of  the  MF  theory  [33]. 
Therefore,  our  calculations  and  the  recent  MC 
studies  of  vacant  sites  requirements  in  adsorption 
of  one  of  the  two  reactants  in  a  bimolecular  sys¬ 
tem  [6]  suggest  that  whenever  a  nonlinear  term 
appears  in  the  ordinary  differential  equations  de¬ 
scribing  the  dynamics  of  the  system,  the  MF  the¬ 
ory  can  be  only  an  approximation.  The  error  in 
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the  approximation  depends  on  the  values  of  the 
parameters,  the  importance  of  migration,  and 
surface  defects. 

MC  simulations  reveal  that,  when  nucleation 
occurs,  the  existence  of  a  Langmuir  process  (reac¬ 
tion  in  the  case  of  equilibrium  transition  and 
desorption  in  the  case  of  kinetic  transition)  de¬ 
stroys  nucleation  and  prevents  the  filling  of  the 
surface  which  causes  poisoning  of  the  surface  re¬ 
action  when  vacant  sites  for  reaction  are  required. 
A  study  of  the  influence  of  migration  on  nuclea¬ 
tion  and  phase  transitions  indicates  that,  if  the 
diffusion  coefficient  is  concentration  independent, 
migration  randomizes  the  distribution  of  atoms 
and  makes  the  MFA  more  accurate.  However,  in 
the  presence  of  attractive  interactions,  migration 
can  facilitate  cluster  formation  and  drive  the  sys¬ 
tem  faster  to  the  stable  steady  state,  making  the 
MFA  worse. 

The  most  severe  limitations  of  the  MF  theory 
concern  its  incapability  to  predict  the  existence  of 
metastable  states  and  the  fluctuations  of  the  sys¬ 
tem  occurring  in  a  small  surface  region.  The  dy¬ 
namic  behavior  of  a  reactive  irreversible  system 
around  a  cusp  point  is  very  important,  because  the 
very  large  fluctuations  can  be  confused  with  the 
oscillations  seen  in  experiments.  The  analysis  of 
dynamics  at  cusp  points  is  important  to  dis¬ 
tinguish  from  the  periodic  patterns  observed  in 
real  experiments.  The  question  of  when  these 
fluctuations  are  related  to  global  periodic  oscilla¬ 
tions  is  investigated  elsewhere  [7], 

In  addition,  no  information  can  be  obtained 
from  the  MF  theory  about  the  influence  of  migra¬ 
tion  on  the  transition  from  the  metastable  state  to 
the  stable  e.  The  influence  of  diffusion  on  kinetic 
transitions  nd  on  Hopf  bifurcation,  as  demon¬ 
strated  elsewhere  [7],  is  totally  missing  from  the 
MF  theory.  Furthermore,  inhomogeneities  in  bi¬ 
nding  energy  with  the  substrate  and  the  existence 
of  point  and  line  defects  which  are  of  great  practi¬ 
cal  importance  in  real  experiments  cannot  be  han¬ 
dled  at  all  by  the  MFA.  Microscopic  models,  such 
as  Monte  Carlo  and  molecular  d  namics  simula¬ 
tions,  in  combination  with  the  new  surface  tech¬ 
niques  could  provide  insight  in  the  dynamics  of 
heterogeneous  processes. 
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